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ABSTRACT 

We report the detection of UCF-1.01, a strong exoplanet candidate with a radius 0.66 ± 0.04 times that of 
Earth (/?©). This sub-Earth-sized planet transits the nearby M-dwarf star GJ 436 with a period of 1.365862 ± 
8x 10~ 6 days. We also report evidence of a 0.65 ± 0.06 R® exoplanet candidate (labeled UCF-1.02) orbiting 
the same star with an undetermined period. Using the Spitzer Space Telescope, we measure the dimming of 
light as the planets pass in front of their parent star to assess their sizes and orbital parameters. If confirmed, 
UCF-1.01 and UCF-1.02 would be called GJ 436c and GJ 436d, respectively, and would be part of the first 
multiple-transiting-planet system outside of the Kepler field. Assuming Earth-like densities of 5.515 g/cm 3 , 
we predict both candidates to have similar masses (^0.28 Earth-masses, Mq, 2.6 Mars-masses) and surface 
gravities of ~0.65 g (where g is the gravity on Earth). UCF-1.01's equilibrium temperature (T ecj , where emitted 
and absorbed radiation balance for an equivalent blackbody) is 860 K, making the planet unlikely to harbor life 
as on Earth. Its weak gravitational field and close proximity to its host star imply that UCF-1.01 is unlikely 
to have retained its original atmosphere; however, a transient atmosphere is possible if recent impacts or tidal 
heating were to supply volatiles to the surface. We also present additional observations of GJ 436b during 
secondary eclipse. The 3.6- \xm light curve shows indications of stellar activity, making a reliable secondary 
eclipse measurement impossible. A second non-detection at 4.5 |xm supports our previous work in which we 
find a methane-deficient and carbon monoxide-rich dayside atmosphere. 
Subject headings: planetary systems — stars: individual: GJ 436 — techniques: photometric 



1. INTRODUCTION 

The search for Earth-sized planets around main-sequence 
stars has progressed expeditiously in the last year. Recent 
discoveries include two Earth-sized planets (0.8 68 and 1.03 
Earth radii, from the Kepler-20 system ( Fressin et alj 
120 121) , two planet candi dates (0.759 and 0.867 R @ ) from the 
KIC 05807616 system dCharpinet et alj|201 lb . and a three- 
planet system (0.78, 0.73, and 0.57 R®) orbiting KOI-961 
ilMuirhead eUdllMLl . 

The search for a second planet in the GJ 436 system be- 
gan shortly after t he transit detection and confirmed eccen tric 
orbit of GJ 436b dGillon et alJl2007t iDeming et alJl2007h . In 
2008, a ^5-M© planet on a 5.2-day orbit was proposed by 
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iRibas et ail (I2008L later retracted) due to three lines of evi- 
dence. First, the lack of detectable GJ 436b transits at the 
time of its 2004 discovery us ing radial-velocity (RV) mea- 
surements (But ler et alj 12004) suggests a change in orbital 
inclination due to a perturb er. Second, given a circulariza- 
tion timescale of ^30 Myr (Demi ng et alJ)2007b and the es- 
timated 6-Gyr age of the system ( Tor res! 12007b . GJ 436b's 
non-circular orbit suggests another planet is pumping up its 
eccentricity. Third, there was evidence of a residual low- 
amplitud e RV signal in a 2:1 mean-motion resonance with 
GJ 436b (IRibas et al.ll2008b . The infe rred planet was discred- 
ited by orbital-dyna mic simulations dBean & Seifahrtl [20081 : 
iDemorv et alj |2009) and the absence of transit timing varia- 
tions (TTVs) with two transit events with the Near Infrared 
Camera and Multi Object Spectrograph camera on the Hub- 
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ble Space Telescope dPont et al J 12009) and ov er a 254-day 
span using ground-based H-band observations dAlonso et alj 
l2O08h. 

iBallard et alJ(l2010bl) 's analysis of 22 days of nearly contin- 
uous observations of GJ 436 during NASA's EPOXI mission 
ruled out transiting exoplanets >2.0 /?® outside GJ 436b's 
2.64-day orbit (out to a period of 8.5 days) and >1.5 in- 
terior to GJ 436b, both with a confidence of 95%. Aided 
by a ^70-hour Spitz er observation of GJ 436 at 8.0 [J.m, 
IBallard et al.l (|2010a) postulated the presence of a 0.75-/?® 
planet with a period of 2.1076 days. However, the predicted 
transit was not detected in an 18-hour follow-up observa- 
tion with Spitzer at 4.5 |_im. The candidate transit signals 
in the EPOXI data w ere likely the result of correlated noise 
(Ball ard et alJl20T0ah . 

In this paper we present Spitzer primary-transit observa- 
tions of UCF-1.01 and UCF-1.02 at 4.5 |om (including an 
independent analysis), a phase curve of GJ 436b at 8.0 \±m 
in which transits of UCF-1.01 are modeled, and a publicly- 
available EPOXI light curve phased to the period of UCF- 
1.01. We also include secondary-eclipse observations of GJ 
436b at 3.6 and 4.5 nm. 

In Section [2] we describe the observations and data analy- 
sis. Section|3]presents Time-series Image Denoising (TIDe, a 
wavelet-based technique used to improve image centers) and 
provides an example analysis using a fake dataset. In Sec- 
tion!?] we discuss the specific steps taken with each of the six 
Spitzer datasets, the details of our independent analysis, and 
transit results from the EPOXI light curve. Section|5]describes 
how we eliminate false positives, our radial-velocity analysis, 
mass constraints on both sub-Earth-sized exoplanets, and or- 
bital and atmospheric constraints on UCF-1.01. Finally, we 
give our conclusions in Section [6] and supply the full set of 
best-fit parameters with uncertainties in the Appendix. 

2. OBSERVATIONS AND DATA ANALYSIS 

2.1. Observations 

We observed GJ 436 at 3.6 and 4.5 |i.m using Spitzer 's In- 
fraRed Array Camera dFazio et al.l 120041) . Including the two 
previously analyzed data sets listed in Table Q] we present six 
Spitzer observations spanning just over three years. 

2.2. POET Pipeline and Modeling 

Our Photometry for Orbits, Eclipses, and Transits (POET) 
pipeline produces systematics-corrected light curves from 
Spitzer Basic Calibrated Data. We flag bad pixels, calcu- 
late image centers fro m a Gaussian fit, and ap ply interpolated 
aperture photometry (Harrin gton et al.l 120071) with a broad 
range of aperture sizes in 0.25-pixel increments. To achieve 
more precise image centers in the 2010 January 28 dataset, we 
utilize TIDe (see Section [3j> . For a more detailed descriptio n 
of POET, see lCampo et allfeoill) and lStevenson et all (1201 11) . 

We model the light curve as follows: 

F( X ,y,t) = F s E(t)R(t)S(t)M(x,y), (1) 

where F(x, y,t) is the measured flux centered at position (x,y) 
on the detector at time t ; F s is the (constant) system flux out- 
side of transit events; E(t) is the primary-transit or secondary- 
eclipse model component; R(t) = 1 -e~ r °' +r ' +r2(t-r^) is the 
time-dependent ramp model component with free parameters 
ro-ry, S(t) = so cos [27r(f -£[)//?] is the phase variation at 8.0 



l_im with free parameters sq and si and p being the fixed pe- 
riod of GJ 436b; and M(x,y) is the Bilinearly-Interpolated 
Subpixel S ensitivity (BLISS) map . We follow the method de- 
scribed by lStevenson et al] (|201 ll) when determining the opti- 
mal bin sizes of the BLISS maps. 

The uniform- sour ce and small-planet equations 
( Mand el & Agoll 120021) describe the secondary-eclipse and 
primary-transit model compon ents. We app l y a non-linear 
stella r limb-darkening model dClaretl [2000; Beaulie uet al.l 
2008) to UCF-1.01 transits with coefficients a\-a^ = 
(0.79660, -1.0250, 0.82228, -0.26800). Spitzer data has well 
documented systematic effects that our Levenberg-Marquardt 
minimizer fits simultaneousl y with the transit/ecl ipse pa- 
rameters. BLISS mapping ( Stev enson et al.l 1201 lb models 
the position-dependent systematics (such as intrapixel vari- 
ability and pixelation) and linear or asymptotically constant 
exponential functions model the time-dependent systematics. 

A Metropolis Random- Walk Markov-chain Monte Carlo 
(MCMC) algorithm assesses the uncertainties ( Cam poet al.l 
l20llb . Each MCMC run begins with a least-squares mini- 
mization, a rescaling of the 5p/teer-supplied uncertainties so 
that the reduced x 2 = L and a second least-squares minimiza- 
tion using the new uncertainties. We test for convergence ev- 
ery 10 5 steps, terminating only when the Gel man & Rubinl 
(1992) diagnostic for all free parameters has dropped to 
within 1 % of unity using all four quarters of the chain. We 
also examine trace and autocorrelation plots of each parame- 
ter to confirm con vergence visually. We estimate the effective 
sample size (ESS. iKass et al.l ll998) and autocorrelation time 
for each free parameter and apply the longest autocorrelation 
time from each event to determine the number of steps be- 
tween independent samples in each MCMC chain. We place 
a prior on UCF-l.Ol's semi-major axis (a/R* = 9.1027^;^) 
by applying its known perio d and GJ 436b's well constrained 
semi-major axis and period (Knutson et al. 201 1) to Kepler's 
third law. Without a prior, the uncertainties for the semi-major 
axis (and any correlated parameters) are larger, but not unsta- 
ble. We also place a flat prior on UCF-1.02's ingress/egress 
time of < 0.1 hours because it is unconstrained by the data. 

3. TIME-SERIES IMAGE DENOISING 

Here we describe an application of wavelets that improves 
image centering, resulting in more precise aperture photome- 
try and better handling of the position-dependent systematics. 
Readers primarily interested in the science results can skip to 
Section|4] 

3.1. Introduction 

Photon noise in short exposures can cause significant shifts 
between the fitted and real stellar centers. With imprecise 
centering over multiple frames, varying amounts of light fall 
within the improperly placed apertures, thus increasing light- 
curve scatter. The sensitivity to precise centering increases 
with smaller aperture sizes, causing a given change in aper- 
ture position to produce larger changes in flux. To improve 
centering, one could sum many exposures, but wavelet filter- 
ing allows the same noise reduction over a short er time span. 
This is important because Stevenso n et al.1 (120101) detect 0.04- 
pixel (0.05 arcsec) pointing variations for IRAC data over ~5 
seconds at > 10cr, which limits the span of an averaging win- 
dow to a few seconds. Our wavelet filter is called Time-series 
Image Denoising (TIDe, pronounced "tidy"). It affects only 
high-frequency components, such as photon noise, without af- 



Two nearby sub-Earth-sized exoplanet candidates in the GJ 436 system 



3 



TABLE 1 
Observation Information 



Observation Date 


Duration 


Frame Time 


Total Frames 


Spitzer 


Wavelength 


Previous 




[minutes] 


[seconds] 




Pipeline 


[nm] 


Publications 1 


2008 July 14 


4,207 


0.4 


588,480 


S18.18.0 


8.0 


K10 


2010 January 28 


1,081 


0.1 


488,960 


S18.18.0 


4.5 


B10 


2010 June 29 


363 


0.4 


49,536 


S18.18.0 


4.5 




2011 January 24 


369 


0.4 


51,712 


S18.18.0 


4.5 




2011 February 1 


369 


0.1 


168,576 


S18.18.0 


3.6 




201 1 July 30 


258 


0.4 


36,160 


S18.18.0 


4.5 





a K10 = parts were published by Knutson et al. I20j7j), B10 = IBallard et aT]<2010al) . 



fecting low-frequency components like transits or eclipses. It 
retains the time resolution of the data. 

In addition to improving aper ture photometry, precise cen- 
tering (see example in Section 13.3b improves our ability to 
model and remove position-dependent systematics accurately, 
for example by reducing the smallest meaningful bin size for 
BLISS mapping. TIDe does correlate the data in time, which 
complicates error analysis and makes it computationally in- 
tense because the correlation depends on the signal and thus 
varies in time. So, we use TIDe-cleaned images only for cen- 
tering (whose uncertainties do not propagate), and perform 
photometry on the unfiltered images. 

As with a windowed (sliding) Fourier transform (WFT), 
wavelets decompose a signal into independent contributions 
at each scale and location (similar to frequency and time) 
within the signal. As an example, the Fourier transform of 
a piece of music can discern the average pitch and timbre of 
all the instruments, but wavelets can identify individual notes 
and the instruments that played them at any given time. The 
wavelet transform of a univariate time series thus has two di- 
mensions, for location and scale. Unlike the WFT, wavelets 
do not suffer from a fixed resolution (or window size), so 
they retain both good temporal resolution for high-frequency 
events and good scale re solution for low-frequency events. 
iTorrence & Compol (11998b provide an accessible introduction 
to wavelets. 

TIDe's improvement in precision and benefits to the light- 
curve fit can vary based on the source brightness, aperture 
size, BLISS map resolution, etc. This method is applicable 
to most photon-noise-limited photometric observations where 
the cadence is significantly shorter than the duration or period 
of the time-varying object of interest. 

3.2. Description of TIDe 

TIDe applies discrete wavelet denoising independently to 
multiple time series, each comprised of the values measured 
in a single pixel as a function of time (i.e., frame number). 
Every pixel is associated with such a time series, and each 
one is denoised independently of adjacent image pixels. The 
transformed data (known as wavelet coefficients) for each 
pixel time series have a location (or time) dimension and 
a scale (or level) dimension. The wavelet coefficients map 
the discrete wavelet to the data at each scale and instant in 
time. The lowest level (or finest scale) of decomposition in 
a wavelet transform describes how the data change on the 
shortest timescales. Assuming that this level is dominated by 
noise, we can eliminate wavelet coefficients with magnitudes 
below a certain threshold (hard thresholding) or merely atten- 
uate them (soft thresholding) to reduce their contribution to 
the overall signal. Adjusting a collection of estimates together 
in this way can be shown to improve the average quality of the 



estimates by introducing a small bias that is more than com- 
pensated for by reduced variance. These techniques can also 
be applied to successively higher levels, but they have less 
impact at longer timescales where the signal dominates over 
the noise. After thresholding, we recombine all of the ad- 
justed wavelets to generate a less-noisy version of the original 
pixel time series. For each frame, an image is re-created from 
the many denoised time series, and centering is performed us- 
ing that image. There is no explicit spatial denoising, but to 
the extent that there are spatial correlations between images 
at different epochs, there is an implicit spatial denoising in 
the processed image that improves center estimation. The ef- 
fectiveness of TIDe is determined by the threshold at which 
wavelet coefficients are zeroed, the type of thresholding tech- 
nique appli ed, and the number of leve l s to which the met hod 
is applied (Donoho & Johnstondll994UChang et aTll2000l) . 

Various wavelet thresholding techniques exist, each with its 
own advantages and disadvantages. Two common methods 
for suppressing noise are hard and soft universal thresholding 
and are defined, respectively, as follows: 

u = yI(\y\>T), (2) 
w = sgn(yXM-r)/(M>7'), (3) 

where y (ui) are the original (denoised) wavelet coefficients 
at a particular level, / is the Indicator function (1 if true, if 
false), and T is some threshold limit. In both instances, if a 
particular wavelet coefficient, Vj, is less than T, then uj i = 0. 
With hard thresholding, the remaining coefficients are unal- 
tered; however, soft thresholding shrinks these coefficients by 
the threshold limit. 

There are many ways to estimate the value of T, including 
VisuS hrink, SURE Shrink, and Bayes Shrink dChang et alj 
2000). With TIDe, we implement the last technique because 
it establishes a thresholding rule that is optimal in terms of 
minimizing the expected RMS error in the denoised time se- 
ries under flexible assumptions for the true time series sig- 
nal (i.e., it minimizes the Bayes Risk for a squared-error 
loss function). Bayes Shrink employs soft thresholding be- 
cause its optimal estimator yields a smaller risk than hard- 
thresholding's estimat or. The optimal thr eshold value is deter- 
mined as follows (see lChang et al.ll2.Q00h . In some instances, 
the noise variance, er 2 , may be known a priori. If this is 
not the case, it is estimated fro m the robust median estima- 
tor dDonoho & Johnstondl 19941) : 

MedianflriQQl) ... 
a= 0.6745 ' (4) 

where Yi(y) represents the wavelet coefficients, y, at the low- 
est level (or finest scale) of decomposition. Next, we esti- 
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mate the variance of Y(y) at a particular level j, assuming zero 
mean, by: 



1 



!=1 



(5) 



where n is the number of wavelet coefficients at that level. 
Our observation model (data = signal + noise) tells us that 
a 2 = a 2 x + a 2 , where a 2 is the signal variance. To account for 
the case where a 2 > a 2 , we calculate a x as follows: 



g x = J max( a 2 - a 2 , 0) 



(6) 



Finally, the optimal threshold at a particular level is deter- 
mined by: 



4. LIGHT-CURVE FITS AND RESULTS 

We present the scaling of the RMS model residuals vs. bin 
size (a test of correlation in time) in Figure |2]for all four 4.5- 
\xm Spitzer observations. Figu re [3] displays our reanalysis of 
GJ 436 data dBallard et al.ll2010al) plus three new Spitzer light 
curves at 4.5 \xm. Two fortuitous detections of UCF-1.01 
appeared during atmospheric ch aracterization observations of 
GJ 436b (Steve nson et al.ll2010t) . Using these data and a ten- 
tative detection at 8.0 [im (see Section |4~6l > to estimate its or- 
bital period, we extrapolated UCF-1.01 transit times forward 
by six months to predict an event (2011 July 30) during the 
next observing window. We supply correlation plots and his- 
tograms in Appendix lAl and the full set of best-fit parameters 
from a 2.4 x 10 6 -iteration joint fit in Appendix IE1 Below, we 
discuss each observation in detail to explain how we arrived 
at the final results. 



In the event that a 2 > a 2 (7) = oo), all of the wavelet coeffi- 
cients are set to zero. 

In this paper, we use the Biorthogonal 5.5 discrete wavelet 
from the PyWavelets package to apply soft thresholding with 
Bayes Shrink to the designated scale levels. 

3.3. TIDe Example 

We generate a series of 1000 test frames, each containing a 
2D Gaussian with a width of 0.7 pixels, a peak flux of 1000, 
and centered at (4.5, 4.5) in a 10 x 10 frame with the lower-left 
corner indexed as (0, 0). We then added a background flux 
offset of 100 and applied Poisson noise to each frame. We 
performed 2D Gaussian centering to derive the blue points 
plotted in Figure Q] For the points in red, we applied TIDe 
to the frames using a Biorthogonal 5.5 discrete wavelet (from 
the PyWavelets package) then recalculated the image centers 
with the same 2D Gaussian centering routine. In each case, 
only the y component of the position is plotted for each frame. 
Using TIDe, the standard deviation in the pointing about the 
true center decreased from 0.019 to 0.01 1 pixels, for a typical 
improvement of ^40%. We see even better results with the 
2010 January 28 data set, where TIDe improved the pointing 
precision by ^70%. 




Bin Size [sec] Bin Size [sec] 



FIG. 2. — Normalized RMS residual flux vs. bin size (in black) for four 4.5- 
u.m light curves. Black vertical lines at each bin size depict lcr uncertainties 
on the RMS residuals (RMS/VlN, where N is the number of bins). The red 
lines show the predicted standard error for Gaussian noise. The dotted and 
dashed lines indicate the scale length of UCF-l.Ol's best-fit transit ingress 
and duration times, respectively. The excess RMS above the red line in the 
top left panel indicates correlated n oise a t timescales near UCF-l.Ol's transit 
duration and is discussed in Section l4~Tl 



0.06 




200 400 600 800 1000 

Frame 



Fig. 1. — Illustrative example of TIDe that compares image centers 
from simulated noisy frames (blue) to their denoised counterparts 
(red). In each case, only the y component of the position is plotted 
relative to the true center. In this typical example, the standard devi- 
ation in the pointing (a measure of precision) decreased from 0.019 
to 0.011 pixels. 



4.1. 2010 January 28 (4.5-\xm Spitzer Observation) 

Spitzer program 541 (Sarah Ballard, P.I.) monitored GJ 436 
continuously for ^18 hours using 0.1 -second exposures. The 
short exposure time allows us to apply TIDe to the lowest 
four levels (L4) of wavelet decomposition (see Section [3}, 
resulting in a maximum affected time resolution of 1 .6 sec- 
onds. In IStevenson et al.l (120101) . we detect Spitzer pointing 
changes on timescales as short as ~5 seconds, longer than 
TIDe's timescale. In calculating image centers with and with- 
out TIDe, we find that the position consistency between con- 
secutive denoised frames improves by more than a factor of 
three, resulting in an RMS of 0.0015 pixels in x and 0.0011 
pixels in y. More precise image centers decrease flux scat- 
ter with smaller aperture sizes and aid the BLISS map in 
modeling the position-dependent systematics. We apply 5x- 
interpolated aperture photometry to the unmodified frames to 
avoid the computationally prohibitive calculation of estimat- 
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FIG. 3. — Four 4.5-|xm Spitzer light curves of GJ 436 with best-fit mod- 
els. The flux values are corrected for systematics, normalized to the system 
brightness, and binned (with \u error bars). Light curves are vertically sep- 
arated for ease of comparison. The single GJ 436b eclipse, four UCF-1.01 
transits, and two UCF-1.02 candidate transits are indicated by the letters b, c, 
and d, respectively. The transits distinguish themselves by their consistency 
in depth and duration. Although UCF-l.Ol's 2010 January 28 best-fit tran- 
sit time is 20 ± 7 minutes earlier than our predicted time (da shed line), the 
parameter's probability distribution is bimodal (see Figure [14) and the other 
peak is only 6 ± 7 minutes early. We quote the median transit time in Tableff] 
to encompass both possibilities. The three remaining UCF-1.01 transit times 
(see TablefS) occur within five minutes of the predicted times and have a typ- 
ical uncertainty of ±3 minutes. The episodic scatter in flux is most likely 
due to stellar activity, which is expected for an M dwarf and seen in many 
observations of this system. Using the non-detection of GJ 436b in the 201 1 
January 24 dataset, we place a 3 it upper limit of 95 ppm on its eclipse depth, 
resulting in a brightness temperature <780 K. This new 4.5-|j.m secondary- 
eclipse observation supports a methane-de ficient and carbon monoxide-rich 
dayside atmosphere i Stevenson et al. 2010). 

ing uncertainties for the denoised frames, which are correlated 
in time. 

Photometry generates consistent transit depths for all tested 
apertures from 1 .25 to 4.50 pixels, but an aperture size of 2.25 
pixels produces the lowest standard deviation of the normal- 
ized residuals (SDNR). We estimate the background flux us- 
ing an annulus from 7 to 15 pixels centered on the star. The 
light curve (see Figure HI exhibits a strong initial increase in 
pixel sensitivity that we do not model (preclip, q = 10,000). 
As with B10, we detect excess flux (possibly due to stellar ac- 
tivity) near BJD 2455225.23 that contributes to the observed 
correlated noise in Figure [2] We note that frames 19,780 
- 19,839, 82,180 - 82,239, 165,380 - 165,439, 419,780 - 
419,839, and 451,780 - 451,839 are shifted horizontally by 
one pixel, so we flag these frames as bad. A probable micro- 
meteor impact at BJD ^2,455,224.976 caused a sudden shift 
in pointing before returning to its original location. Simulta- 
neously, the background scatter increased by ^50% and re- 



mained elevated until the end of the observation. We apply 
a BLISS map bin size (x, y) of 0.007 x 0.005 pixels and set 
the minimum number of points per bin to six to disregard the 
observed excursion. 




-0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 



BJD - 2455225 

Fig. 4. — Full light curve from 2010 January 28 Spitzer observa- 
tion with transits of UCF-1.01 (right) and UCF-1.02 (left). The flux 
values are corrected for systematics, normalized to the system bright- 
ness, and binned (with la error bars). The solid line depicts the best- 
fit model. Excluding the excess flux near 2455225.23 in the model 
fit does not significantly alter the best-fit solution. 



4.2. 2070 June 29 (4.5-yun Spitzer Observation) 

Our Spitzer program 60003 (Joseph Harrington, P.I.) mon- 
itored GJ 436 for six hours using 0.4-second exposures. The 
mean image center is located at pixel (15, 25), near the top 
of the 32x32 array, thus restricting aperture sizes to <5.50 
pixels. Using a background sky annulus from 10 to 30 pixels, 
we find that the lowest SDNR occurs with a 5 x -interpolated 
aperture 5.00 pixels in radius. The BLISS map uses bins of 
size 0.006 x 0.009 pixels and with at least four points per 
bin. We test image centers generated from L3 TIDe (3.2- 
second maximum time resolution) but find no improvement 
in the SDNR. This is likely due to the smaller improvement 
in image centers and significantly larger aperture size, rela- 
tive to the 2010 January 28 dataset. The final analysis did not 
use TIDe. For this dataset, the telescope pointing does not 
stabilize until midway through the transit of UCF-1.02 (see 
Figure 0. As a result, the position-dependent systematic is 
poorly constrained during and prior to the transit. This may 
be the source of UCF-1.02's variable transit duration, which 
decreases with smaller aperture sizes. More observations are 
necessary to confirm its parameters. 

4.3. 2077 January 24 (4.5-\xm Spitzer Observation) 

Spitzer program 60003 performed a second six-hour ob- 
servation of GJ 436 with 0.4-second exposures. We find 
that 10 x -interpolated aperture photometry outperforms 5x- 
interpolated and minimizes SDNR with an aperture size of 
5.25 pixels and a background sky annulus from 10 to 30 pix- 
els. We flag 54 frames (28,426 - 28,479) as bad due to a 
one-pixel horizontal shift, as observed previously in a dataset 
above. We clip the first 6,000 observations due to a strong 
increase in flux, possibly due to stellar activity (see Figure[3]i. 
Near 2455585.771, we observe a sudden shift in the telescope 
pointing that, again, correlates with an increase in background 
noise. To remove this excursion from our models, the BLISS 
map uses bins with eight or more points and a size of 0.016 
x 0.008 pixels. As with the previous dataset and for the same 
reasons, TIDe centers have little effect on the resulting photo- 
metric light curve. 

4.4. 2077 February 1 (3.6-\un Spitzer Observation) 



6 



Stevenson et al. 



15.05 

X 

c 15.00 
c 

.3 14.95 

£ 14.90 
~Q) 

X 14.85 

Q_ 

14.80 
26.0 

>. 

c 25.8 
c 

.2 25.6 
"ui 

£ 25.4 
>< 25.2 

Q_ 

25.0 










UCF-1.01 
H UCF-1.02 








to/ft 


***** 



0.60 0.65 0.70 

BJD - 2455225 



Binned Data 
Best Fit 
GJ 436b 



I Hi 



0.54 0.56 0.58 0.60 0.62 

Orbital Phase (2.64-day period) 

Fig. 6. — Spitzer 3.6-(im light curve from 2011 February 1 with 
GJ 436b's time of secondary eclipse shaded in gray. The flux val- 
ues are corrected for systematics, normalized to the system bright- 
ness, and binned (with la error bars). The solid line depicts the 
best-fit baseline model. Stellar activity prohibits us from fitting the 
eclipse and measuring its depth. However, by visually comparing the 
binned points within the shaded region to those immediately outside, 
the data appear t o be consistent with a relatively deep eclipse depth 
(Stevensoh'et al. 2010). 




14.86 14.88 14.90 14.92 14.94 14.96 14.98 15.00 

Pixel Position in x 

Fig. 5. — Image centers vs. time (upper 2 panels) and pointing his- 
togram (lower panel, number of image centers within a given bin) 
for the 2010 June 29 dataset. The times during transit are shaded in 
gray. Initial telescope drift hampers our ability to effectively model 
position-dependent systematics during and prior to the UCF-1.02 
transit. 

Our Spitzer program 60003 also observed GJ 436 at 3.6 \im 
with 0.1 -second exposures. We apply 5 x -interpolated aper- 
ture photometry with an aperture size of 2.75 pixels and a 
background sky annulus from 7 to 15 pixels. We clip the 
first 10,000 frames due to a steep ramp and frames 70,000 
- 125,000 due to suspected stellar activity (see Figure|6]l. Be- 
cause GJ 436b's time of secondary eclipse occurs during the 
period of increased stellar activity, we do not fit the eclipse or 
calculate uncertainties. 

4.5. 2077 July 30 (4.5-\xm Spitzer Observation) 

Spitzer monitored GJ 436 for 4.3 hours using 0.4-second 
exposures (program 70084, Joseph Harrington, P.I.). Pho- 
tometry generates consistent transit depths for apertures be- 
tween 1.75 and 6.00 pixels. The final run applies 10 x- 
interpolated, 5.00-pixel aperture photometry with a back- 
ground sky annulus from 10 to 30 pixels. During these ob- 
servations, the telescope pointing experiences two deviations, 
at BJD 2,455,772.766 and 2,455,772.870. The background 
variance increases with the first event but slightly decreases 
with the second event. BLISS mapping utilizes a bin size of 
0.012 x 0.006 pixels with a minimum of six points per bin 
to exclude points from either excursion. Again, TIDe centers 



have little effect on the resulting photometric light curve. 

4.6. 2008 July 14 (8.0-\xm Spitzer Observation) 

Spitzer program 50056 (Heather Knutson, P.I.) observed GJ 
436 for -70 hours from 2008 July 12 to 2008 July 15. At 
the best aperture size of 3.75 pixels (and a background sky 
annulus from 7 to 15 pixels), we find that the light curve ex- 
hibits a meas urable position-depend ent systematic, identified 
as pixelation (IStevenson et al.ll201 ll) . The BLISS map fits and 
removes pixelation (see Figure |7]l using a bin size of 0.022 
x 0.022 pixels and at least four points per bin. We model 
the initial time-dependent ramp using an asymptotically con- 
stant exponential function after clipping the first 3,000 frames. 
A sinusoidal function with a linear correction fits the phase 
variation of GJ 436b (see Figure [S). We set a prior on the 
inclination and semi-major axis of UCF-1.01 using the best- 
fit results from the 4.5-(xm joint fit. The UCF-1.01 transit at 
BJD 2,454,662.328 i s the same candidate transit reported by 
iBallard et al.l (|20 1 0a) using a ~2 . 1 -day period estimated from 
EPOXI observations. Their Figure 5 incorrectly reports the 
BJD. We used the timing of this transit to successfully predict 
the 2001 July 30 transit. The best-fit radius ratio from both 
transit events in this light curve is 0.010 ± 0.003, which is 
consistent with the best-fit result using the four 4.5- [J.m light 
curves. 

4.7. Independent Analysis 

We sought an independent analysis to confirm our results. 
Nikole Lewis analyzed each of the 4.5- |am datasets without 
knowing the times or depths of the transits. In addition to 
using her own photometric pipeline, she applied a new pixel- 
mapping rout i ne that shares a heritage with the method from 
Ballar d et al.l d2010al) . This new pixel-mapping method was 
developed to recover the relative flux variations as a func- 
tion of orbital phase from the Spitzer 3.6-ixm and 4.5-iim 
full orbit light curves of HD 189733b, HD 209458b, HAT- 
P-2b, and HAT-P-7b (PLKnutson; PID 60021). Similar to 
the BLISS method, the pixel-mapping technique developed 
by Lewis uses nearest neighbors to calculate flux as a func- 
tion of position on the detector, but in her method the dis- 
tances are weighted according to a Gaussian distribution. In 
addition to stellar centroid positions, Lewis makes use of the 
"noise pixel" parameter given in frame headers to determine 
the nearest-neighbors to a given data point (Lewis et al., in 
prep.). This routine improves on Ballard's method by calcu- 



Two nearby sub-Earth-sized exoplanet candidates in the GJ 436 system 



7 



15.0 




Pixel Position in x 




Pixel Position in x 



Fig. 7. — BLISS map (top) and pointing histogram (bottom) for the 
2008 July 14 dataset. Pixelation, a position-dependent systematic, is 
depicted by the colors in the BLISS map, where redder (bluer) col- 
ors indicate more (less) flux within the aperture. Peaks repeat every 
0.2 pixels because we applied 5 x -interpolated aperture photometry. 
Smaller interpolation factors result in larger spacing between peaks 
but also a stronger systematic between peaks. The horizontal and 
vertical black lines depict pixel boundaries. 

lating the pixel map at each iteration without being computa- 
tionally prohibitive. 

Pixel mapping is essential to detecting the weak transit sig- 
nals in these data. For example, the 2010 January 28 dataset 
requires an accurate pixel mapping routine, at minimum, to 
detect UCF-1.01 and benefits from more precise image cen- 
ters with TIDe to more clearly distinguish UCF-1.02. We 
have found that without a pixel-mapping routine, one can- 
not reproduce all of the observed transits. Lewis uncovered 
transits of UCF-1.01 in the 2010 June 29, 2011 January 24, 
and 2011 July 30 datasets with ease and, once informed of 
the additional planet, identified both UCF-1.02 transits and 
the remaining UCF-1.01 transit (see Figure |9). Her final tran- 
sit times, depths, and durations for both planets are all within 
1.5a of our best-fit results. 

4.8. EP OXl Observation 

NASA's EPOXI mission obser ved GJ 436 nearly co ntinu- 
ously during 2008 May 5-29 dBallard et al.l l2010b). The 
light-curve data are available from EPOXI's archive. After 
masking the transits of GJ 436b, we divide the light curve 
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Fig. 8. — Spitzer light curve of GJ436 at 8.0 (im. The flux values are 
corrected for systematics, normalized to the median system bright- 
ness, and binned (with la error bars). The solid blue line depicts the 
best-fit model. The light curve contains two eclipses (first and last 
events) and one transit (center event) of GJ 436b and two transits 
(second and fourth events) of UCF-1.01. The two UCF-1.01 tran- 
sit depths have a combined significance of 2a, which is insufficient 
to claim a detection, but we used the timing of the latter to predict 
the 2001 July 30 transit. The difference in brightness temperatures 
between GJ 436b's dayside and nightside causes the observed sinu- 
soidal variation in the light curve. The peak-to-peak flux difference 
is 200 ± 50 ppm (4a significance). This corresponds to a brightness 
temperature difference of 110 ± 60 K, which favors a relatively ef- 
ficient dayside-to-nightside energy redistribution. The peak flux is 
shifted by 0.7 ± 4.6 hours prior to secondary eclipse. 

by the median flux value, phase it according to the best-fit 
UCF-1.01 period (see Tabled, and bin the results. Figure [TOl 
illustrates a visible decrease in the observed flux at the correct 
phase that is consistent with the transit depth and duration of 
UCF-1.01 derived from the Spitzer data. The quality of the 
light curve is such that the data neither prove nor disprove the 
existence of UCF- 1.01. 

5. DISCUSSION 

Without continuous monitoring of GJ 436 for two consecu- 
tive transits at the most photometrically-precise wavelengths 
(3.6 and 4.5 (a.m), we isolate the true period from integer mul- 
tiples or whole number fractions by other means. Integer mul- 
tiples (i.e., 2, 3, 4...) of the orbital period (see Table|2| cannot 
account for all of the observed transits; whole number frac- 
tions (i.e., 1/2, 1/3, 1/4...) are eliminated by investigating 
the bevy of available GJ 436 Spitzer data at predicted tran- 
sit times. We find evidence against periods of ^0.6829 and 
-0.4553 days by non-detections of UCF-1.01 in a 2008 Jan- 
uary 30 observation at 3.6 pin and a 2008 June 1 1 observation 
at 8.0 \xm, respectively. The single UCF-1.01 detection in the 
2010 January 28, 18-hour observation dismisses even shorter 
periods. 

5.1. Eliminating False Positives 

GJ 436's large proper motion (across the sky) enables us to 
eliminate astrophysical false positives that could mimic the 
observed periodic decrease in flux. Over our 1.5-year ob- 
servational baseline of 4.5- jam detections, the system moves 
~1.8", equivalent to 1.5 pixels in Spitzer 's InfraRed Array 
Camera. With aperture sizes as small as 1.25 pixels for the 
first (2010 January 28) and last (2011 July 30) observations, 
we find that the transit signals from UCF-1.01 are clearly dis- 
tinguished against the background noise. This limits the lo- 
cation of a potential background sourc e (such as an eclipsing- 
binary star system, Torre s et al]|201 11) to the overlapping re- 
gion within both apertures. Using obser vations from the Very 
Large Telescope (V LT, see Figure 

03] iRousset et all [2003; 
iLenzen et al.1 120031) an d Canada France Hawaii Telescope 
(CFHT, see Figure [L2l iRigaut et al.1 1 19981) with adaptive op- 
tics imaging instruments at two different epochs, we elimi- 
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Fig. 9. — Four 4.5-nm Spitzer light curves of GJ 436 with best- 
fit models from an independent analysis by Lewis. She corrects 
flux measurements for intrapixel sensitivity variations using a pixel- 
mapping technique and for the presence of the well documented 
Spitzer time-dependent systematic (ramp). A fixed- width symmet- 
ric Gaussian fits centroid positions in the region near the bright- 
est pixel in each subarray frame. The best photometric aperture 
size is 2.25 pixels for the 2010 January 28 dataset and 5.0 pixels 
for the other datasets. The non-linear limb-dark ening coefficients 
for GJ 436 are those from Knu tson et al.l d2011l) . After the loca- 
tion of the transit(s) in each dataset were identified individually, 
Lewis performed a simultaneous fit between all four datasets using 
a Levenberg-Marquardt minimization scheme. A MCMC algorithm 
determined the uncertainty in the fit parameters as well as identified 
other possible solutions. The goal of this analysis was to confirm 
the presence and shape of transit(s) in each dataset. Improvements 
to treatment of the systematics in these observations is possible, but 
they are unlikely to significantly change the estimated planetary pa- 
rameters. 




Fig. 10.— EPOXI light curve phased to the period of UCF-1.01 
using the best-fit period and nearest ephemeris time (2008 July 14 
dataset). Blue circles represent the binned EPOXI data with lrj un- 
certainties. The red cross depicts the duration and depth (with a io 
uncertainty) of UCF-1.01's transit. The EPOXI data are consistent 
with a UCF-1.01 transit. 

nate background stars up to 12.7 and 9.3 magnitudes fainter, 
respectively, than GJ 436 at a confidence of 5a. 



TABLE 2 

Transit model best-fit values and other parameters 



Parameters 


UCF-1.01 


UCF-1.02 


Rp/R* 


0.0138 ± 0.0009° 


0.0136 ±0.0012 




85.17^" 


- 


a/R* 


9.10 ± 0.07" 


- 


Impact Parameter 


u.; i_q 15 


- 


Transit depth [ppm] 


190 ± 25 


186 ± 30" 


Duration P4.i1 hr] 


76*0- 15 


1 04+0.26n 


Ingress/Egress [hr] 


025+0.002 


<0.1° 


Transit Times [MJDtto] 6 


5225.090^0°^ 


5225.026 ± 0.003" 




5376.7078+o ; oo''* fl 


5376.568^™ fl 




5585.6889 + o ; oo20" 






5772.8069^;0°<» fl 




Mean Period [Days] 


1.365862 ± 8'xlfr 6 




Ephemeris [MJD7-0B]* 


5772.8086 ± 0.0016 




Radius [Rgj] 


0.66 ± 0.04 


0.65 ± 0.06 


Mass [M 9 ] d 


0.28 ± 0.06 


0.27 ± 0.07 



Fitted values. 



b MJD = BJD - 2,450,000. 

' We choose the median value because the distribution is bimodal. 
d Assuming an Earth-like density of 5.515 g/cm 3 . 




FIG. 11. — Very Large Telescope H-band observation o n 2007 March 20 
using the NAOS-CONICA instrument with adaptive optics I Montagnier et al. 
2012). We search for faint background systems by blocking the light from GJ 
436 using a 0.7" Lyot coronagraphic mask. The dark green lines are mask 
support wires. The "+" symbols indicate the position of the GJ 436 system for 
this observation and at each transit epoch of UCF-1.01. Red circles indicate 
the minimum photometric aperture size ( 1 .25 pixels) for which transit signals 
from the first and last confirmed events may still be clearly distinguished 
against the background noise. If a background system were the source of the 
transit-like events, it must put light in the overlapping region. To produce 
the observed transit depth, the hypothetical system must be no more than 
9.3 magnitudes fainter than GJ 436, assuming a total eclipse of one of the 
objects. We eliminate objects brighter than AH = 12.7 relative to GJ 436 with 
a confidence of 5cr. There are also no objects listed in any catalog within this 
region. 

5.2. Instability Hypothesis 

In this section we consider an alternate hypothesis to that 
of detecting two sub-Earth-sized exoplanets, that the observed 
variations are the result of instrumental or stellar instabilities. 
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FIG. 12. — Canada France Hawaii Telescope (CFHT) K-band observation 
obtained 2000 April 19 using the adaptive optics bonnette (PUEO). The "+" 
symbols indicate the position of the GJ 436 system for this observation and at 
each transit epoch of UCF-1.01. Red circles indicate the minimum photomet- 
ric aperture size (1.25 pixels) for which transit signals from the first and last 
confirmed events may still be clearly distinguished against the background 
noise. We eliminate background objects within the overlapping region with 
AK = 9.3 at a 5cr confidence limit. 

We begin by calculating the probability of observing UCF- 
1.01 and UCF-1.02 by chance. Then, we quantify the occur- 
rence rate of transit-like instabilities and estimate the proba- 
bility that these instabilities are periodic. Finally, we compare 
our model fits to the null hypothesis, which is expressed by 
a model that does not contain planet parameters, to see if the 
additional free parameters are justified. 

In search of transit signals in the GJ 436 system, we ex- 
amined 11 light curves at 3.6 and 4.5 [am (not counting the 
201 1 July 30 dataset in which we predicted the transit). Both 
channels have the photometric stability necessary to detect 
GJ 436c. Of the 71.3 hours of data, there are eight transit 
or eclipse events of GJ 436b, each lasting ^-4 -hour in dura- 
tion. Since we cannot distinguish overlapping transits, we 
have 63.3 hours of usable data with an average light-curve 
duration of ^5.75 hours. Given that the period of UCF-1.01 
is 1.3659 days, the probability that a transit will occur in a 
typical event is 17.5%. Using the binomial distribution, we 
calculate a 30.1% chance of observing three or more transits 
of UCF-1.01 in the 11 available light curves. Recall that our 
fourth transit event of UCF-1.01 was predicted, rather than 
occurring by chance, so it does not enter into the calculation. 

We repeat the above calculation for UCF- 1 .02 but must first 
estimate its orbital period by considering the transit duration 
ratio between itself and GJ 436b. We find that the durations 
are nearly identical; however, both planets are unlikely to oc- 
cupy the same orbit. So, we perform two sets of calculations: 
one for each side of the 1-sigma uncertainty in UCF-1.02's 
transit duration. Using a period of 5.563 days, the probability 
of observing two or more UCF-1.02 transits is 7.9%. Using 
a period of 1.785 days, the probability of observing two or 
more transits increases to 44.6%. The combined probability 
of observing both planets ranges from 2.3% to 13.4%. 



We compare these results to the alternate hypothesis, 
namely that the observed variations are the result of instru- 
mental or stellar instabilities. To begin, we analyzed nearly 
120 hours ofGJ 1214data at 4.5 |im (Sptizer program 70049). 
This M dwarf is similar to GJ 436 and should exhibit similar 
levels of activity (stellar instabilities). If the instabilities are 
instrumental, then it should not matter which star we analyze 
unless the instabilities are flux-dependent (GJ 436 is almost 3 
magnitudes brighter than GJ 1214 in the infrared). From our 
GJ 1214 light-curve results, we identify two transit-like events 
based on their depths (>200 ppm) and durations (~1 hour). 
Assuming these events are not the result of planet transits, for 
any given hour of 4.5-|im observations, there is a 1.7% prob- 
ability of having an instability event. We then apply the bino- 
mial distribution to determine that the probability of detecting 
five or more instabilities in 63.3 hours of data is 0.42%. Re- 
call that we do not count the 201 1 July 30 dataset or use times 
during GJ 436b transits/eclipses. If we assume that the insta- 
bilities only appear at 4.5 microns, the probability of detect- 
ing five or more instabilities in 44.7 hours of data decreases 
to 0.088%. 

Since we cannot find a physical mechanism for reducing the 
stellar flux in a transit-like way with a repeatable period, any 
observed instabilities must be random events. We consider 
the probability of detecting four random instability events that 
happen to coincide with a given period (in this case, 1.3659 
days). The first two instability events establish the "period" 
under consideration. As calculated above, the third and fourth 
instability events each have a 1.7% probability of occurring 
within 30 minutes of the established period (total time win- 
dow is 1 hour). Their combined probability is 0.029%. 

We conclude that the single-planet hypothesis is 72 times 
more likely than the most favorable instrumental/stellar- 
instability scenario and 1040 times more likely than detect- 
ing four random instability events that happen to coincide 
with a given period. The two-planet hypothesis is 5.5 - 32 
times more likely than the most favorable instrumental/stellar- 
instability scenario and 79 - 460 times more likely than de- 
tecting four random instability events that happen to coincide 
with a given period. 

Finally, we test the strength of our two sub-Earth-sized exo- 
planet candidates by comparing various fits to the null hypoth- 
esis, which asserts that there are no new planets. Recall that 
a lower ABIC value indicates that the additional free parame- 
ters are warranted in the model fit. Relative to the null hypoth- 
esis, ABIC decreases by 11.4 when we include UCF-1.01's 
transit parameters in a joint model fit. Alternatively, if we 
add only UCF-1.02's transit parameters then ABIC increases 
by 36.2. Including both planets' transit parameters in a joint 
model fit results in an increase in ABIC of 14.5, relative to the 
null hypothesis. Thus, the BIC favors a model that includes 
UCF-1.01 but disfavors models that include UCF-1.02. This 
result is directly related to the number of observed transits 
for each planet candidate and indicates that we need to obtain 
more than two transit observations of UCF-1.02 to increase 
the detection significance and surpass the BIC's penalty for 
additional free parameters. We conclude that the available 
data support UCF-1.01 as a strong exoplanet candidate and 
signify that UCF-1.02 is a weak exoplanet candidate. 

5.3. Radial-Velocity Constraints 

The 3.6-meter ESQ telescope at La Silla Observatory 
dMavor et al.1 120031; iPepe et alj 120041) utilized the HARPS 
spectrograph with the settings described by iBonfils et alj 
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( 12011b to obtain 171 observations of GJ 436 at 550 nm. 
Xavier Bonfils (personal communication) provided us with 
the extracted, unpublished RV measurements so that we 
could attempt to constrain the mass of UCF-1.01. We re- 
tained 159 data points (12 were removed due to the Rossiter- 
McLaughlin effect). To these data, we ad ded 41 GJ 436b 
primary transit times (Knutson et ail 1201 ll and references 
therein), 14 GJ 436b secon dary eclipse times dStevenson et alj 
l2010t iKnutson et alj|201 ll), and an 8.0- \im photometric light 
curve from lDeming et aiT(l2007l) . The light curve (retrieved 
from the Infrared Processing and Analysis Center, IPAC) is 
binned into 445 points and normalized to remove the time- 
dependent ramp. 

We apply a two-planet model with the transit ephemeris for 
the second planet fixed to the best-fit value listed in Table [2] 
and the eccentricity fixed to 0. The fit utilizes the empirical 
stellar density calibration of lEnoch et alJ (120101) to determine 
the stellar mass, in addition to other system parameters, grant- 
in g this fit a much broa der scope than the modeling described 
bv lCampo et al.l (1201 ll) . We employ a Levenberg-Marquardt 
minimizer to find the best-fit parameters to our model and a 
Markov-chain Monte Carlo routine with 10 6 iterations to esti- 
mate uncertainties. We express our \ 2 function as follows: 



Vj-Vj 

0V,i 



-E 



! 2 



+ E 

k 



Ok~ Ok 



0~ ,k 



■E 
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(8) 



where v, t, o, and p represent the HARPS radial velocities, 
primary-transit times, secondary-eclipse times, and photomet- 
ric data, respectively. The over-lined quantities indicate com- 
puted values and a represents the uncertainty for each mea- 
surement. We adjust for transit-eclipse light travel times and 
for leap seconds in this fit. Using the above data with an es- 
timated stellar jitter of 1.7 m/s, we do not detect the signal of 
the second planet but cannot repudiate its existence. The 3a 
upper limit of the semi-amplitude is 0.6 m/s, corresponding 
to an upper limit of 0.6 M®, which is larger than our mass 
constraints using the density arguments below. 

5.4. Mass Constraints 

Unable to effectively constrain the mass of UCF- 1.01 using 
RV data, we consider a range of possible bulk densities for 
a terrestrial-sized planet (see Figure [T3l>. Given a mean bulk 
density between 3 and 8 g cm -3 , we limit the mass of UCF- 
1.01 to be 0.15 - 0.40 M ffi and estimate the surface gravity 
to be 0.36 - 0.94 g. We place similar limits on the mass and 
surface gravity of UCF- 1.02. Assuming an Earth-like density 
of 5.515 g cm" 3 , we estimate masses of 0.28 and 0.27 M ffi for 
UCF-1.01 and UCF-1.02, respectively, which correspond to 
surface gravities of ^0.65 times that on Earth. 

5.5. Orbital Constraints 

UCF-1.01 may exhibit TTVs due to gravitational interac- 
tions with GJ 436b in a near-2:l orbital resonance or with 
UCF-1.02, which has an unknown orbit. This may explain 
why UCF-1.01's transit time is 20 minutes early in the 2010 
January 28 data set; however, the parameter's probability dis- 
tribution is bimodal (see Figure[T4li and the other peak is only 
6 ± 7 minutes early. The three remaining UCF-1.01 transit 
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FIG. 13. — Mass and surface-gravity constraints on UCF-1.01 (solid lines) 
and UCF-1.02 (dashed lines). Bold lines depict the best-fit values and thin 
lines depict the upper and lower lcr uncertainties. Exoplanet KOI-961.03 and 
solar-system planets are included for reference. 

times occur within five minutes of their predicted times and 
have a typical uncertainty of ±3 minutes. More precise ob- 
servations could establish whether these deviations are TTVs. 

Using the known orbital parameters of GJ 436b and UCF- 
1.01, we performed orbital-stabili ty simulati ons using the 
Mercury numerical integrator ( Chambers 1999). Assuming an 
Earth-like density of 5.5 g cm" 3 , the predicted mass of UCF- 
1.01 is 0.28 Mm. We supplied the code with initial starting 
conditions, listed in Table [3] based on transit times from the 
2011 January 24 dataset. Our results indicate that the orbits 
are stable out to at least 100 Myr. The best-fit line shows a 
change in semi-major axis of 5.3e-06 au/Gyr. The osculat- 
ing UCF-1.01 orbital parameters exhibit a periodic trend ev- 
ery ~35 years wherein the eccentricity varies between and 
0.21, the peak-to-trough inclination amplitude is 3.2°, and 
TTVs vary from ±200 to ±3 minutes. A ~40-day periodic 
trend is also evident but with smaller variations in the oscu- 
lating orbital parameters. Due to UCF-1.01's relatively small 
mass, variations in GJ 436b's orbital parameters over the 35- 
year timespan are below Spitzer's sensitivity limits. Next- 
generation facilities may be able to constrain UCF-l.Ol's or- 
bital parameters through improved RV measurements or by 
measuring its time of secondary eclipse. 

5.6. Atmospheric Constraints 

UCF-1.01 is unlikely to have retained any original atmo- 
sphere due to its weak gravitational field, close proximity to 
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TABLE 3 

Initial orbital parameters for GJ 436b and UCF-1.01 



Parameter 


GJ 436b 


UCF-1.01 


Semi-major Axis 


0.0287 au 


0.0185 au 


Eccentricity 


0.1371 





Inclination 


86.699° 


85.1° 


Argument of Periapsis 


351.0° 


0.0° 


Longitude of the Ascending Node 


0.0° 


0.0° 


Mean Anomaly 


282.6° 


90° 



its ho st star, and estimated 6-Gyr age of the system dTorresI 
120071) . The planet receives a substantial soft x-ray and ex- 
treme ultr aviolet (XUV) flux; we estimate 700 - 900 er g 
cm -2 s _1 dSanz-Forcada et al.ll20"TTt lEhrenreich et all 1201 11) , 
or — 1,000 times the present XUV flux received by the Earth. 
Such an intense XUV flux leads to a v ery hot the rmos phere 
and subsequent hydrodynamic escape (Tian 2009). Shortly 
after formation, outgassing from an Earth-like, silicate-rich 
mantle could have pr oduced an initial wate r- vapor-rich atmo- 
sphere for UCF-1.01 (Schaefe r et al.1201 11) . However, the wa- 
ter vapor would readily have been photolyzed by ultraviolet 
radiation at high altitudes, leading to a hydrogen-dominated 
thermosphere that like ly extended to the planet's Roche dis- 
tance of -25,000 km dErkaev et alj|2007l) . given the planet's 
low gravity. In this situatio n, the mass-loss ra te for energy- 
limited hydrodynamic flow dErkaev et alJl2007l) implies a hy- 
drogen loss rate of about 8 x 10 10 g s" 1 (assuming an XUV 
heating efficiency of 1), or 1.4 times the planet's mass lost 
in 1 Gyr. This indicates that hydrogen was lost from UCF- 
1 .01 's atmosphere very early in its history. Some heavy ele- 
ments would have been entrained in the hydrodynamic flow, 
but the early atmosphere would have become increasingly ox- 
idized as hydrogen was lost. Carbon dioxide could then have 
dominated at some later point in the atmosphere's history, but 
even a C02-rich atmosphere would be unstable. Scaling from 
hydrodynamic models (Tian 2009), we estimate that carbon 
would be lost from a C02-rich atmosphere at — 1 x 10 8 g 
s , or 1% of the planet's mass over its lifetime - an amount 
likely greater than the planet's initial inventory of CO2. At- 
mospheres dominated by molecula r nitrogen o r oxygen would 
be lost on even shorter timescales (Tian 2009). 

UCF-1.01 could support a transient, present-day atmo- 
sphere if recent impacts were to deliver volatiles rather than 
preferentially erode any atmosphere, or if tidal heating were 
to supply volatiles from the crust/mantle. The latter sce- 
nario is particularly attractive if a recycling mechanism ex- 
ists for any heavy atmospheric constituents (e.g., volcanic 
emission of sulfur dioxide, followed by photolysis to sulfur 
and oxygen atoms, dayside-to-nightside transport, condensa- 
tion, and subsequent melting and re-vaporization of sulfur de- 
posits). In this speculative scenario, UCF-1.01 could resem- 
ble a hot Io that has lost its lighter and more volatile elements. 



Any transient atmosphere will likely have a low surface pres- 
sure and be highly extended, which could fill the Roche lobe 
and/or produce a tail. Transit observations at ultraviolet wave- 
lengths could confirm or rule out such an extended atmo- 
sphere, and one might search particularly at wavelengths in 
which atomic and ionized sulfur and oxygen would be ex- 
pected to absorb. Given that volcanically supplied sodium 
and potassium might be transient atmospheric constituents, 
visible-wavelength transit observations might also prove use- 
ful. 

6. CONCLUSIONS 

In this paper, we announced the detection of UCF-1.01 
and UCF-1.02, two sub-Earth-sized transiting exoplanet can- 
didates orbiting the nearby M dwarf GJ 436. Their detec- 
tions were possible with BLISS mapping and Time-series Im- 
age Denoising (TIDe), the latter of which is a novel wavelet- 
based technique that decreases high-frequency noise in short- 
cadence, time-series images to improve image centering pre- 
cision. We presented four transits of UCF-1.01 and two tran- 
sits of UCF- 1 .02 at 4.5 \im, an independent analysis that con- 
firms our best-fit results within 1.5<r, an 8.0-Lim phase curve 
of GJ 436b that includes transits of UCF-1.01, and EPOXI 
data that are consistent with the presence of a sub-Earth-sized 
exoplanet. To definitively establish UCF-1.01 as a planet (to 
be called GJ 436c), we require only a few hours of addi- 
tional observations, preferably from another telescope or at 
least at a different wavelength. Establishing UCF-1.02 as a 
planet (to be called GJ 436d) would likely require an extended 
observing campaign to constrain its period then successfully 
predict a transit. Fi nally, we confirmed the GJ 436b 4.5- |am 
results presented by Stev enson et al.l (120101) through an addi- 
tional non-detection during secondary eclipse; however, we 
were unable to confirm the strong eclipse depth at 3.6 [am due 
to stellar activity. The current data still support a methane- 
deficient and carbon monoxide-rich dayside atmosphere. 
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Fig. 14. — Correlation plots and histograms for the 18-hour 2010 January 28 Spitzer observation containing transits of UCF-1.01 and UCF- 
1.02. We plot every 4000 th step in the MCMC chain to decorrelate parameter values. UCF-l.Ol's distribution of mid-transit times (midpoints) 
is bimodal, so we favor the median value over the best-fit value (see Table [2j- UCF-1.02's ingress/egress times are unconstrained from our 
model fit. 
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Fig. 15. — Correlation plots and histograms for the 2010 June 29 Spitzer observation containing transits of UCF-1.01 and UCF-1.02. We plot 
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Fig. 16. — Correlation plots and histograms for the 201 1 January 24 Spitzer observation containing a transit of UCF-1.01 and an eclipse of GJ 
436b. We plot every 4000 th step in the MCMC chain to decorrelate parameter values. 
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B. BEST-FIT PARAMETERS 



TABLE 4 

Best joint-fit light-curve parameters 



Parameter 


2010 January 28 


2010 June 29 


2011 January 24 


2011 July 30 


2008 July 14 


Wavelength ( ixm) 


4.5 


4.5 


4.5 


4.5 


8.0 


Array Position (x, pix) 


14.69 


14.94 


14.63 


14.70 


14.54 


Array Position (_v, pix) 


14.92 


25.31 


15.20 


14.98 


14.52 


Position Consistency' 1 (<5 V , pix) 


0.0015 


0.0025 


0.0045 


0.0041 


0.0055 


Position Consistency' 1 (S y , pix) 


0.0011 


0.0039 


0.0028 


0.0024 


0.0052 


Aperture Size (pix) 


2.25 


5.00 


5.25 


5.00 


3.75 


Inner Sky Annulus (pix) 


7.0 


10.0 


10.0 


10.0 


7.0 


Outer Sky Annulus (pix) 


15.0 


30.0 


30.0 


30.0 


15.0 


System Flux, F s (|xJy) 


841090 ± 100 


871540 ±30 


819510 ±30 


825590 ± 15 


315195 ±7 


GJ 436b Tr. Midpt." (MJDtdb) 


- 


- 


- 


- 


4661.50365 ± 0.00012 


GJ 436b «„/«„ 


- 


- 


- 


- 


0.0830 ± 0.0006 


GJ 436b cosi 


- 


- 


- 


- 


0.066 ±0.002 


GJ 436b a/R, 


- 


- 


- 


- 


13.0 ±0.3 


GJ 436b Eel. Midpt." (MJDtdb) 


- 


- 


- 


- 


4660.417 ±0.003 


GJ 436b Eel. Midpt." (MJDtdb) 


- 


- 


5585.7747 


- 


4663.053 ± 0.003 


GJ 436b Eel. Duration fa-i, hrs) 


- 


- 


1.00 


- 


1 .02 ± 0. 1 3 


GJ 436b Eclipse Depth (ppm) 


- 


- 


18 ±28 


- 


500 ±60 


GJ 436b T b (K) 


- 


- 


540 ± 80 


- 


700 ± 30 


GJ 436b Amplitude, .50 (ppm) 


- 


- 


- 


- 


100 ± 40 


GJ 436b Offset 6 , si (MJDtdb) 


- 


- 


- 


- 


4660.39 ±0.19 


UCF-1.01 Midpt." (MJDtdb) 


5225.090!g;™ 


5376.7078!g;S 


5585.6889^ooofs 


5772.8069^™' 


4662.328 ±0.013 


UCF-1.01 Rp/R* 


0.0138 ± 0.0009 


0.0138 ± 0.0009 


0.0138 ± 0.0009 


0.0138 ± 0.0009 


0.010 ±0.003 


UCF-1.01 cos/ 


0.084_^ 


084+ 0003 


0.084tJg 


0.084*"" 


0.084 ±0.003 


UCF-1.01 a/R* 


9. 10 ±0.07 


9.10 ±0.07 


9. 10 ±0.07 


9. 10 ±0.07 


9. 10 ±0.06 


UCF-1.02 Midpt." (MJDtdb) 


5225.026 ± 0.003 


5376.568!U;J!o7 


- 


- 


- 


UCF-1.02 Transit Depth (ppm) 


186 ±30 


186 ±30 


- 


- 


- 


UCF-1.02 Duration ((4-1, hrs) 




1 04+ - 22 


- 


- 


- 


UCF-1.02 Ingress (f 2 -i, hrs) 


0.06 ±0.03 


0.06 ±0.03 


- 


- 


- 


UCF-1.02 Egress (ft- 3 , hrs) 


0.06 ±0.03 


0.06 ±0.03 


- 


- 


- 


Ramp, ro 


7.0 ±2.0 





22 ± 12 


44 ±11 


18.1 ± 1.6 


Ramp, n 


-8.4 ± 0.6 





6±8 


25 ±8 


-0.6 ± 0.5 


Ramp, r2 


-0.0008 ± 0.0003 


0.0020 ± 0.0003 








-0.00012 ± 0.00006 


Ramp, ri 





0.5 








1.5 


TIDe 


Yes 


No 


No 


No 


No 


BLISS Map [M(x,y)] 


Yes 


Yes 


Yes 


Yes 


Yes 


Effective Sample Size (ESS) 


341 


702 


398 


675 


415 


Minimum # of Points Per Bin 


6 


4 


8 


6 


4 


Total Frames 


488960 


49536 


51712 


36160 


588480 


Rejected Frames (%) 


0.178 


0.527 


0.673 


0.465 


0.393 


Frames Used 


477106 


48777 


44728 


35172 


583049 


Free Parameters 


6 


10 


5 


4 


18 


AIC Value 


605808 


605808 


605807 


605808 


583067 


BIC Value 


606091 


606091 


606090 


606091 


583270 


SDNR 


0.00535600 


0.00253643 


0.00257029 


0.00258144 


0.00508140 


Uncertainty Scaling Factor 


0.31734 


0.17676 


1.06515 


0.98102 


1.04102 


Photon-Limited S/N (%) 


84.3 


82.2 


84.0 


83.2 


84.3 



"RMS frame-to-frame position difference. 
b MJD = BJD - 2,450,000. 

c We exclude frames during instrument/telescope settling, for insufficient points at a given knot, and for bad pixels in the photometry aperture. 



